Pair contact process with diffusion of pairs 
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Abstract 

The pair contact process (PCP) is a nonequilibrium stochastic model which, hke the basic contact 
process (CP), exhibits a phase transition to an absorbing state. The two models belong to the 
directed percolation (DP) universality class, despite the fact that the PCP possesses infinitely 
many absorbing configurations whereas the CP has but one. The critical behavior of the PCP 
with hopping by particles (PCPD) is as yet unclear. Here we study a version of the PCP in which 
nearest-neighbor particle pairs can hop but individual particles cannot. Using quasistationary 
simulations for three values of the diffusion probability {D = 0.1, 0.5 and 0.9), we find convincing 
evidence of DP-like critical behavior. 
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I. INTRODUCTION 



The exploration of universality classes associated with nonequilibrium phase transitions 
continues to attract much interest l|-[3| . In the case of systems exhibiting a phase transitions 
;o an absorbing state, the generic universality class is that of directed percolation (DP) 
4, sl- For example, both the basic contact process (CP) jsl and the pair contact process 
(PCP) 0, Q belong to this class, despite the fact that the former has a unique absorbing 
configuration while the latter possesses infinitely many. Allowing individual particles to hop 
in the PCP, one obtains the so-called pair contact process with diffusion (PCPD). In this 
case there are only two absorbing states: the empty lattice, and the state of a single particle 
hopping. 

A model with the basic features of the PCPD was first proposed by Grassberger in 1982. 
More recent interest in this problem was stimulated by the work of Howard and Tauber 9| 
who discussed a generalized version, with a L ang evin description involving complex noise. 
On the basis of numerical results, Carlon et al. 10|, suggested that the critical behavior of the 



PCPD would fall in the parity-conserving (PC) class. Subsequent works [ll|, ll2| suggested 
a different kind of critical behavior, possibly masked by huge corrections to scaling. (For a 
comprehensive review of analyses of the PCPD up to 2004, see the article by Henkel and 



Hinrichsen 



13|. 



Given the controversies surrounding the PCPD, it seems worth checking whether allowing 
hopping by pairs {only) changes the critical behavior of the PCP. One expects such a model 
to remain in the DP universality class, given that the critical behavior of the CP is robust 
to the inclusion of nearest-neighbor hopping [l^. Note as well that pairwise diffusion does 
not change the set of absorbing configurations in the PCP. In the present work, we study 
a variant of the PCP in which nearest neighbor pairs of particles may hop together, while 
isolated particles cannot hop. Using quasistationary (QS) simulations, we determine several 
critical parameters for three values of the diffusion probability: D = 0.1,0.5 and 0.9. The 
balance of this paper is organized as follows. In Sec. II we define the model and detail our 
simulation method. In Sec. Ill we present our results; Sec. IV is devoted to discussion and 
conclusions. 



2 



II. MODEL AND SIMULATION METHOD 

The CP O] is one of the simplest and most studied models belonging to the DP 
universality class. In the CP, each site z of a lattice is either occupied [cilt) = 1] or vacant 
[o"j(t) = 0]. Transitions from cTj = 1 to occur at a rate of unity, independent of the 
neighboring sites. The reverse transition can only occur if at least one neighbor is occupied: 
the transition from (jj = to 1 occurs at rate Am, where m is the fraction of nearest neighbors 
of site i that are occupied; thus the state = for all i is absorbing. A is the only control 
parameter in the basic CP; the order parameter p is the fraction of occupied sites. 



The PCP 



8| is defined on a rf-dimensional lattice, with each site again either occupied 
or vacant. Only pairs of particles occupying nearest-neighbor sites exhibit activity: each 
such pair has a rate 1 — p of mutual annihilation, and a rate p to create a new particle 



3, 



at a randomly chosen site neighboring the pair, if this site is vacant. In the PCPD 
in addition to the creation and annihilation processes present in the PCP, each particle 
attempts to hop, at rate D, to a randomly chosen nearest-neighbor (NN) site; the move is 
accepted if the target site is vacant. Several variants of the model, diffe ring in how each 



process (creation, annihilation or diffusion) is selected, have been studied |13 |. 

The model studied here, pair contact process with diffusion of pairs (PCPDP), again 
features the pair-mediated annihilation and creation processes of the PCP. In addition, a 
fraction D of all events are hopping attempts by a NN particle pair. The model is defined 
on a ring of L sites. The transition rates associated with each pair (AA) are: 

dlAA AAdS, AAdl^dlAA, at rate L'/2 

AAdi, USAA^AAA, at rate p(l - D)/2 (1) 

AA^m, at rate (1 - £)). 

Thus hopping ceases, along with all other activity, in the absence of pairs. 

The most obvious definition of the order parameter in the PCPDP is the pair density 
Pp, that is, the number of nearest- neighbor occupied pairs divided by the total number of 
nearest-neighbor pairs on the lattice (on the ring, the latter is equal to the number of sites). 
Since the number of particles can only change in the presence of pairs, one might expect the 
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excess particle density, Appart = Ppart — Pnnrt,c, to scale in a similar manner. (Note that the 
critical particle density Ppart,c is nonzero |7|.) We confirm that Appart and pp exhibit similar 
scaling properties. 

We sample the quasistationary (QS) distribution of the process, i.e., conditioned on 
survival, using a simulation method that yields quasistationary properties directly |l5 |. 
This is done by maintaining, and gradually updating, a set of configurations visited during 
the evolution; when a transition to the absorbing state is imminent the systems is instead 
placed in one of the saved configurations. Otherwise the evolution is identical to that of 
a conventional simulation. The set of saved configurations is updated by replacing with a 
small probability, Prepi at each time step, one of the saved configurations with the current 
one. 

We perform extensive simulations of the one-dimensional PCPDP, using systems of L = 
10,20,40,80,160,320,640 sites, using the QS simulation method. Each realization of the 
process is initialized with all sites occupied, and runs for 5 x 10^ to 10^ time steps (longer 
runs for larger systems). Our results are calculated over samples sizes 10-20 realizations. 
The number of saved configurations ranges from 100 to 1000 (larger numbers for smaller 
systems). Values of p^ep range from 10^^ to 5 x 10"^ (smaller values for larger systems). 
During the relaxation phase we use a prep ten times larger, to eliminate the influence of 
the initial configuration. Following this relaxation phase, we accumulate histograms of the 
time during which the system has exactly 1, 2, n, ... pairs, and similarly for particles. 
Using the histograms we calculate the densities of pairs and of particles, the moment ratio 
m = (Pp) / {Pp)"^ and the reduced fourth cumulant = K^/ K^, where K2 = var(pp) and, 

= {pt) - MpI){Pp) - HpI? + 12(pJ)(Pp)^ - Q{Pp)\ (2) 

The QS lifetime r is taken as the mean time between attempts to visit an absorbing config- 
uration. 



III. SIMULATION RESULTS 



For each diffusion rate studied, the first step in the analysis is to determine the critical 
annihilation probability Pc{D). Experience with absorbing-state phase transitions leads us 
to expect the following scaling properties at the critical point: and 
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m(L) —J- rric, a universal critical value 16|. We obtain a preliminary estimate of pc from the 
crossing of the moment ratios m for system sizes L = 10 and 20, and then concentrate our 
efforts in simulations near this value. (For D = 0.5, for example, we focus on the interval 
p G [0.85,0.90], see Fig. [I]). 
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FIG. 1: (Color online) Moment ratio m versus creation parameter p for D = 0.5, L = 10, 20. Inset: 
Reduced fourth cumulant versus p for L = 160, D = 0.5. 

For each pair of consecutive lattice sizes, L and 2L, we determine the crossing values px 
and Tfix- An independent estimate of the critical point is afforded by the reduced fourth 
cumulant, q^, which takes a pronounced minimum at a value Pq{L) that converges to the 
critical value; see Fig. [H inset. In Fig. |2]we plot p^ and Pq versus 1/L, for D = 0.5; both 
sets of values converge quickly to a limit which we estimate, using quadratic extrapolation 
(versus l/L^), as p^ = 0.8900(2). 

The finite-size scaling relations for pp, r, and m, cited above, are satisfied to good precision 
for p = 0.8900; the data for the order parameter are shown in the inset of Fig. |2l Including 
results for L = 640, however, we find small but significant curvatures in the graphs of In pp 
and Inr versus InL, indicating that Pc is slightly larger than 0.8900. (The curvatures are: 
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FIG. 2: (Color online) Values px and pq associated, respectively, with moment ratio crossings 
(lower) and minima in the reduced fourth cumulant (upper), versus 1/L, for D = 0.5. Inset: Order 
parameter pp versus system size, for D = 0.5 and p = 0.89. 

-0.0019(4) in the case of Pp and -0.009(3) for r.) We perform additional simulations for 
L = 640 at p = 0.8903 and 0.8906, and then calculate the curvatures of the graphs of In pp 
and Inr versus InL for diverse values of p in the interval [0.89, 0.8903], using polynomial 
interpolation of the simulation data as necessary, to estimate the quantities of interest at 
intermediate points. We observe no significant curvature for p in the interval [0.89004, 
0.89010], leading to our final estimate, p^ = 0.89007(3), for D = 0.5. The associated critical 
parameters are (3/i^± = 0.252(2), z = 1.573(10), and m = 1.1758(24). (Note that the 
principal source of uncertainty in these estimates is the uncertainty in pc itself.) 

The data for the particle density Ppart scales in the same manner as the pair density: 
at the critical point, we find, to good precision, ppart — Ppart,c + AL~'^^'^^, using the same 
value for as found in the analysis of the pair density. (For D = 0.5, for example, 

the limiting particle density is Ppart,c = 0.2321(2).) The moment ratio associated with the 
particle density, nipart = {{Ppart- PpaTt,c)'^) / {Ppart - Ppart,c)^ , behaves similarly to the moment 
ratio associated with pairs; for D = 0.5 we find mpart,c = 1.177(5). 
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FIG. 3: (Color online) Derivatives \dm/dp\ (lower), dlnr/dp (upper), and dlnpp/dp (middle) 
versus L for D = 0.5. 

The simulation data also permit a direct estimate of the exponent z/^, albeit with some- 
what limited precision. Finite-size scaling implies that the derivatives \dm/dp\, dlnr/dp 
and dlnpp/dp, evaluated at the critical point, follow \dx/dp\ oc L^^'^^ (here x stands for any 
of the quantities mentioned). We estimate the derivatives via least-squares linear fits to the 
data on an interval that includes pc- (The intervals are small enough that the graphs show 
no significant curvature.) Power-law dependence of the derivatives on system size is verified 
in Fig. 131 Linear fits to the data for the three largest sizes, for m. In pp, and In r yield 
l/u^ = 0.900(4), 0.940(26) e 0.891(7), respectively, leading to the estimate u± = 1.10(2). 
The simulations and analyses described above were repeated for diffusion rates D = 0.1 
and 0.9, yielding the results shown in Table I. For D = 0.9, finite-size corrections are much 
stronger than for the smaller diffusion rates, and we found it necessary to extend the study 
to larger system sizes to obtain reliable results. The values listed in Table I are based on 
studies using L = 320, 640, 1280, and 2560, in this case. 
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D 


yc 




z 


V 1 


' ' ''p 


0.1 


0.91720(1) 


0.252(1) 


1.584(11) 


1.11(4) 


1.173(1) 


0.5 


0.89007(3) 


0.252(2) 


1.573(10) 


1.10(2) 


1.1758(2) 


0.9 


0.83470(5) 


0.253(1) 


1.573(5) 


1.12(3) 


1.178(2) 


DP 




0.25208(5) 


1.5807(1) 


1.096854(4) 


1.1736(1) 



TABLE I: Critical exponent values for the PCPDP and DP; figures in parentheses denote uncer- 



tainties. DP values from Refs. [la ] and [1 

IV. CONCLUSIONS 

We study a version of the pair contact process in which nearest-neighbor particle pairs 
(but not isolated particles) diffuse. The results of quasistationary simulations, for (S/ux, 
z, and the moment ratio m, confirm to good precision that the scaling properties of 
this model coincide with those of directed percolation. This conclusion is what one would 
expect on the basis of universality: diffusion of pairs changes nothing fundamental in the 
propagation of activity or in the space of absorbing configurations of the original pair contact 
process. The situation is quite different from that of the PCP with diffusion of particles: 
the critical behavior of the PCPD is still not fully understood, and is evidently subject to 
much stronger corrections to scaling than is the model studied here (PCPDP). Thus our 
conclusion that the PCPDP belongs to the directed percolation universality class does not 
imply that the same holds (or does not hold) for the PCPD. 
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